Partial Weyl Law for Billiards 
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For two-dimensional quantum billiards we derive the partial Weyl law, i.e. the average density of 
states, for a subset of eigenstates concentrating on an invariant region T of phase space. The leading 
term is proportional to the area of the billiard times the phase-space fraction of F. The boundary 
term is proportional to the fraction of the boundary where parallel trajectories belong to V. Our 
result is numerically confirmed for the mushroom billiard and the generic cosine billiard, where we 
count the number of chaotic and regular states, and for the elliptical billiard, where we consider 
rotating and oscillating states. 
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Waves restricted to finite resonators in two or three 
dimensions have found abundant applications in con- 
temporary physics. Examples include electromagnetic 
and acoustic resonators, microdisk lasers, atomic matter 
waves in optical billiards, and quantum dots [IHBj]. The 
average density of states is an essential observable of a 
resonator and dictates many physical properties. Inves- 
tigations of this quantity go back to Lord Rayleigh more 
than a hundred years ago [6| and have been a continuing 
topic of interest ever since [7Hl4j . Today the fundamental 
result for the average density of states 
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is known as Weyl's law, as he gave the first proof of 
the leading term jjjj. Equation (pj is formulated for the 
case of a two-dimensional quantum billiard with area A, 
boundary length L, and Dirichlet boundary conditions 
(in units h — 2m — 1). Analogous results for three di- 
mensions, other types of waves and boundary conditions 
are available [FJ], and for open systems a fractal Weyl 
law was proposed [TEj . 

The first term of eq. ([IJ depends on the area A of the 
billiard only. As for any quantum system it is obtained 
by counting the number of Planck cells in the phase space 
available at energy E. The second term was already con- 
jectured by Weyl @. It is specific for billiards or res- 
onators and depends on the length L of the boundary. 
As the wave function must vanish on the boundary, a 
layer is depleted which has a length L and a width of the 
order of the wave length A ~ 1/ \f~E~. Scmiclassically, the 
second term can be interpreted as a contribution from 
closed trajectories which are reflected perpendicularly at 
the boundary [ToL fl6l IFm] . Higher order corrections in 
eq. CD arise, e.g., due to corners and curvature effects 



Generic billiards, which are studied for electromag- 
netic, acoustic, and matter waves [lHj], have a phase 
space containing several dynamically separated domains, 
such as regular and chaotic regions, see fig. [He). The 



spectrum consists of sub-spectra with eigenfunctions 
mainly concentrating on one of these invariant regions 
Fi, according to the semiclassical eigenfunction hypoth- 
esis It is a fundamental question to know the 
corresponding partial average density of states dr^E), 
where d(E) = ^^dr^E). This is essential for study- 
ing the spectral statistics [2l[ of such sub-spectra. The 
partial density of states is also compulsory for the deter- 
mination of transition rates with Fermi's golden rule [22j . 
Furthermore, it is required when an external coupling to 
the system is not uniform in phase space, e.g., the total 
internal reflection in optical resonators [23, [24| or tilted 
leads attached to a quantum dot [25j . 

In this paper we derive the partial Weyl law, 
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for a subset of eigenstates corresponding to an invariant 
region T of phase space, using the Wigner-Weyl trans- 
formation of the Green function. As expected, for the 
area Ar the ratio Ay /A is the fraction of phase space 
occupied by T, see eq. ([F4|. For the length Lr we find 
that the ratio Lr/L is the fraction of the billiard bound- 
ary where parallel trajectories belong to T, see eq. (|15l) . 
This is unexpected as semiclassically the boundary term 
in eq. (QJ originates from trajectories perpendicular to 
the boundary [lfj, [l6|, |l7[ . We confirm the result eq. ([2]) 
numerically for the mushroom billiard and the generic 
cosine billiard, where we predict the number of regular 
and chaotic states, and for the elliptical billiard, where 
we consider rotating and oscillating states. 

Before we derive eq. ^ we exemplify its application 
for the desymmetrized mushroom billiard (26j shown in 
fig. QJc). It is characterized by the radius R of the 
quarter circular cap, the stem width a, and the stem 
height I. Classically one has two distinct phase-space 
regions visualized in fig. [He). All trajectories which 
are located only in the cap of the mushroom are reg- 
ular, while those entering the stem are chaotic [26l j . 
For the chaotic region eq. (fT4")) yields A c h — ttR 2 /4 — 
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[R 2 arcsin(a/i?) + aV 'R 2 - a 2 ]/ 2 [22J. For the length L ch 
we have to consider those parts of the billiard boundary 
9f2 for which parallel trajectories belong to the chaotic 
region (including the marginally stable bouncing-ball or- 
bits). This is the case for the straight boundaries, such 
that L c h = 2i? + 21 is the length of the chaotic bound- 
ary of the mushroom billiard. A TCg and L rcg follow from 
A lcg — A — ^4 c h and L rcg = L — L c \ x = irR/2, respectively. 
(For the full mushroom billiard one finds L rcg = ttR 
and L c h = 2i? + 21.) In order to verify this predic- 
tion of the partial density of states, eq. @, we numer- 
ically solve the time-independent Schrodinger equation, 
— A-0;(q) = Eiipi(q), for the dcsymmctrizcd mushroom 
billiard with Dirichlet boundary condition (ipi(q) = 0, 
q e dft), R = 1, I = 1, and a = 0.5. We calcu- 
late the first 6024 eigenstates using the improved 
method of particular solutions [27]. They can be clas- 
sified as mainly regular or mainly chaotic, depending 
on the phase-space region on which they concentrate 
(fig. [He)). There are several methods for this classifi- 
cation which give similar results. Here we determine the 
regular fraction w xes of an eigenstate Vz of the mush- 
room by its projection on a basis of the regular region. 
For this basis we use the eigenstates ip™ n of a quarter 
circle of radius R = 1 with energy E mn and angular mo- 
mentum m > a\JE mn . These basis states are given by 
$£"( r i V>) = N mn J m (j mn r) sin(m^), where m = 2,4, . . . 
is the angular quantum number, n = 1, 2, . . . is the radial 
quantum number, J m is the mth Bessel function of the 



first kind, j mn is the nth root of J m , E„ 



Jrnn ) an d 

N mn — y/8/ir/ Jm-i (jmn ) is a normalization constant. 
The projection of ipi onto these basis states leads to the 
regular fraction w l rcg = J2 m , n K^"l^>| 2 0(m-a^/E mn ) 
„l 



with < w l reg < 1. The chaotic fraction is then given 
1 - 



and w l ch the densities 



by w l ch = 1 — «4eg- From w\ 
of the regular and the chaotic states can be computed, 
dr(E) — wj, S(E — Ei). For comparison with numer- 
ics we use the more convenient spectral staircase function 
Nr (E) — drj dr (rj) . It has a step of size w l T at eigenen- 
ergy L/. From eq. © one finds 



N r (E) = ^E-^VE. 

47T 47T 
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Figure [TJa) shows the regular and the chaotic spectral 
staircase for the mushroom billiard. We find excellent 
agreement with our prediction, eq. © (smooth solid 
lines). In fig. [TJb) we demonstrate that the boundary 
contribution of eq. ([3]) is in agreement with the differ- 
ence of the numerical data and the first term of eq. ((3|). 
Using the semiclassical interpretation of the boundary 
term in eq. ([1]) one would naively expect that Lr/L is 
the fraction of the boundary where perpendicular trajec- 
tories belong to T. However, this does not reproduce the 
data (dotted lines) . The numerical fluctuations arise due 
to oscillatory contributions to the density of states which 
are not considered here. We have confirmed that under 
variation of the width a of the stem of the mushroom the 
prediction eq. © agrees with numerics. 
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FIG. 1. (Color online) (a) Spectral staircase Nr(E) for reg- 
ular and chaotic eigenstates of the desymmetrized mushroom 
billiard with R — 1, I = 1, and a = 0.5. We compare numeri- 
cal data with the first term (dashed lines) and both terms of 
eq. © (smooth solid lines). The inset shows the regular and 
chaotic parts of the boundary, (b) Same data after subtract- 
ing the area term, Ar(-E) = Nr{E) — ArE/ (4-7r), compared to 
the second term of eq. © (smooth solid lines) shown over a 
larger energy range. Determining Lr in eq. © from perpen- 
dicular instead of parallel trajectories gives an incorrect result 
(dotted lines), (c) Phase space at the circular boundary with 
regular (lines) and chaotic (dots) regions and illustrations of 
trajectories and eigenfunctions. 



Now we turn to the derivation of eq. ©. The most gen- 
eral method to represent quantum states in phase space 
is their Wigner distribution. Other options, such as the 
Husimi distribution or numerical methods for specific ge- 
ometries as in the example above, yield similar results 
and can be obtained from averages over the Wigner dis- 
tribution. In the semiclassical limit, according to the 
semiclassical eigenfunction hypothesis, the weights w l r for 
an invariant region T of non-zero measure are either zero 
or one and independent of the projection method. 

Based on the Wigner distribution of an eigenstate ipi, 



W,(q,p) = ^ y d ; 



re 



2ipr 



^*(q + r)<Mq-r), (4) 



we define the phase-space resolved density of states for a 
phase-space point (q, p) 



d w (£) = ^^-£ 1 )^(q,p) 



(5) 



For an arbitrary region T of phase space the density of 
states is given by the integral 



^(L)=yd 2 qd 2 pd q , p (L). 



(6) 
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Using 



ImG(x,x' )J E) = -7r£;vr(x')^W^-^) (7) 
1=1 

eq. ([U can be rewritten as 

dq, P (£) = -^3 J d 2 re 2ipr lmG(q-r,q + r,E), (8) 

where G(x,x',i?) is the energy dependent Green func- 
tion. For billiards it satisfies 

(A + k 2 )G(x,x',£;) = <J(x-x'), x,x'gO, (9) 

where k = yE. On dfl it satisfies the boundary condi- 
tion of the billiard and is zero outside f2. Close to the 
boundary the curvature can be neglected and the Green 
function of a half-plane is an appropriate approximation. 
It is then given in the upper half-plane by 

G(x,x',£) « ±{H+( K \ X x'|) - H+(k\x x'|)]. (10) 

Here Hq denotes the Hankel function of the first kind 
and x' is the mirror image of x' . 

Instead of requiring that the Green function vanishes 
in the lower half-plane, it is more appropriate to continue 
the Green function antisymmetrically across the bound- 
ary, i.e. eq. (fTU|) is extended to the full plane. In this 
way the corresponding Wigner function is adapted to the 
Dirichlet boundary condition and yields a more faithful 
momentum distribution [28]. Far from the boundary this 
modification has no effect and reduces to the standard 
definition. 

Using polar coordinates for the momentum p = (p, /3) 
and Cartesian coordinates for the position q = {q\\,qx), 
measured parallel and perpendicular to the boundary <9f2, 
we evaluate in eq. © the integrals over r and in eq. © 
the integral over p, 
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]2„„2ipr 



Jo(2/s|r|) - Jo fa^Tqlj] , (11) 



where Jo is the Bessel function of the first kind. For 
the r-intcgration of the first Bessel function we use 
polar coordinates r = (r, tp), where ip is the angle 
with respect to p, § dip e 2ipr cos f = 27rJ (2pr), and 
/ °° drr J a (2pr)J {2nr) — S(p — n)/Ap. For the re- 
integration of the second Bessel function in eq. (fTTj) 
we use Cartesian coordinates r = (r\\,r±) leading to 
pr = pr\\ cos(j3 — + pr± sin(/3 — where (3\\(q) is 
the angle of the tangent to the boundary at q € <90. In- 
tegration over rj_ gives 7n5(sin[/3 — P\\])/p and integration 
over p leads to 7nJ(rii)/2. Finally we obtain 
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1 
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The (5-function selects trajectories with momentum par- 
allel to the boundary. According to eq. (j6|) we have 



ME) 



d 2 q 



(13) 



where xr(q, 0) is the characteristic function of the phase- 
space region T. It is one when the trajectory running at 
angle j3 through the point q belongs to T, zero if this is 
not the case, and 1/2 on the boundary of T. We now 
evaluate eq. (Ti~3"|) in the semiclassical limit k — > oo, where 
Jo(2nq±) — > 5(q±)/K. This gives the final result eq. ([2]) 
with 



A r = f d 2 q±- r d^ X r(q,/3), 
Jn 27T Jo 



Lr = 



ds xr(q(s),/3||(s)). 



(14) 
(15) 



6(wa[/3 - )3 ll ])J (2Kq ± ). (12) 



Here, s is the arc length along the boundary, q(s) is 
the corresponding point on the boundary, and /3ii (s) = 
/3ii(q(s)) is the angle of the tangent to the boundary at 
that point. We thus find that Lr/L is the fraction of 
the billiard boundary where parallel trajectories belong 
to T. Strictly speaking, xr(q(s), /3y (s)) in eq. (fT5t is ob- 
tained as a limit from trajectories starting at Q|| = s with 
—> and q± —> 0. This infinitesimal neighborhood 
has to be considered when the parallel trajectory with 
— /3|| and q± = cannot be assigned to one of the re- 
gions r. Note, that Ar is proportional to the volume of 
r, while there is no such simple relation for Lr- For the 
special case where T is the entire phase space we have 
Xv = 1 leading to Ay = A and Lr = L, such that eq. ([2]) 
reduces to eq. (JTJ) - 

We now give two explanations for the boundary term 
in eq. ([5]) and its relation to trajectories of T that are 
parallel to the boundary: (i) If T is the entire phase 
space eqs. © and §8§ lead to dq(E) = — Im G(q, q, E)/ir. 
The semiclassical contributions to this Green function arc 
given by closed trajectories which start and end at q [l[ . 
According to the second term in eq. (fTU|) they have a 
length 2q± corresponding to a perpendicular reflection 
at the boundary. This leads to the term Jo(2«g^)/(47r), 
which agrees with the second term in eq. (1121) when inte- 
grated over all j3. If T is a region in phase space, one needs 
the phase-space resolved density of states dq :P (E). Equa- 
tion © shows that the direction of the momentum p is 
not related to the direction of the trajectories which semi- 
classically contribute to the Green function. Evaluating 
eq. (fTTj) leads to rii = 0, such that still perpendicularly re- 
flected trajectories of length 2q± give the boundary term. 
However, this contribution arises only if the direction f3 
of the momentum p is parallel to the boundary. For the 
other directions the contribution cancels due to the phase 
factor exp(2ipr) in eq. (fTT|). (ii) An intuitive explanation 
of the boundary term can be given in terms of a plane 
wave, exp(ip||g||) s'm(p±q±), with Dirichlet boundary con- 
dition at q± — 0. Here the sine-term suppresses waves 
with p± — and thus reduces the density of states for 
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these waves, which have a momentum p parallel to the 
boundary. Semiclassically such waves correspond to tra- 
jectories parallel to the boundary, in agreement with our 
result for the boundary term, eq. (fFoj) . 

We now consider the desymmetrized cosine billiard in 
order to demonstrate that the partial Weyl law can be ap- 
plied to systems with a generic mixed phase-space struc- 
ture. The cosine billiard is characterized by the height h 
and length a of the rectangular part as well as the height 
h c of the upper cosine boundary, see fig. [UJc). For the 
chosen parameters the phase-space of the cosine billiard 
consists of one large regular region surrounded by chaotic 
motion and a four-island resonance chain, see fig. [2fc), 
which also shows a magnification of the generic hierar- 
chical regular-to-chaotic transition region. In order to 
apply the partial Weyl law we have to define a region 
r in phase space and determine the corresponding area 
Ar and length Ly- According to eq. © this gives a pre- 
diction for the number of eigenstates concentrating on 
r. In principle any invariant region can be used. We 
choose the red-shaded region in fig. [2Jc) , which contains 
most of the central regular island. The area Ar is de- 
termined from eq. (|14l) by numerical integration. For 
the length Lr we have to consider those parts of the bil- 
liard boundary dfl for which parallel trajectories belong 
to T. This holds for the left vertical boundary of length 
Ly = h + h c , see the inset in fig. Ufa). We stress, that 
including the hierarchical region or parts thereof in the 



definition of T affects Ar, but not the boundary term Lr, 
as the orbits parallel to the boundary do not belong to 
the hierarchical transition region. Numerically we cal- 
culate the first 1853 eigenstates of the desymmetrized 
cosine billiard with h — 0.8, a = 1.3, and h r = 0.24 using 



the improved method of particular solutions 27]. For the 
/th eigenstate we integrate its Poincare-Husimi distribu- 
tion [30j over the region corresponding to T which gives 
the weight w l r . These weights determine the spectral 
staircase function Nr(E). We find excellent agreement 
with our prediction, eq. (J3|), see fig. [^a). In fig.^b) we 
demonstrate that the boundary contribution of eq. ([3]) is 
in agreement with the difference of the numerical data 
and the first term of eq. ©, apart from a constant offset 
due to higher order terms neglected in eq. ([3]). 

As an interesting application, where a part of the 
boundary contributes to two invariant regions of phase 
space , we now consider the desymmetrized elliptical bil- 
liard [31| , shown in fig. |3Jc) . It is characterized by the 
lengths r a and rj, of the two half-axes, with r a > rj,, and 
the focus / = yjr\ — r% . The phase space of the elliptical 
billiard consists of two separated regions of rotating and 
oscillating motion as visualized in fig. [He). The quan- 
tum eigenstates can be classified accordingly as mainly 
rotating or oscillating, N(E) = N Iot (E) + N osc (E). Us- 
ing eq. (j3j we predict the number of rotating and oscil- 
lating states up to energy E. The areas A mt and A osc 
are determined from eq. (I14[) by numerical integration. 
For the length L rot we have to consider those parts of 
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FIG. 2. (Color online) (a) Spectral staircase Nr{E) for eigen- 
states concentrated in an invariant region F of the desym- 
metrized cosine billiard with a = 1.3, h = 0.8, and h c = 0.24. 
We compare numerical data with the first term (dashed line) 
and both terms of eq. © (smooth solid line) . The inset shows 
the part of the boundary corresponding to T. (b) Same data 
after subtracting the area term, Ar(-E) = Nr(E) — ArE/(4n), 
compared to the second term of eq. © (smooth solid line) 
shown over a larger energy range, (c) Phase space at the lower 
horizontal boundary with invariant region V (red shaded), re- 
maining phase space (gray lines and dots), magnification of 
the hierarchical regular-to-chaotic transition region, and illus- 
tration of a regular trajectory. 




FIG. 3. (Color online) (a) Spectral staircase Nr(E) for rotat- 
ing and oscillating eigenstates of the desymmetrized elliptical 
billiard with r a — 1 and n = 0.7. We compare numerical data 
with the first term (dashed lines) and both terms of eq. © 
(smooth solid lines). The inset shows the rotating and oscil- 
lating parts of the boundary, (b) Same data after subtracting 
the area term, Ar(E) = Nr(E) — ArE/(4n), compared to 
the second term of eq. © (smooth solid lines) shown over a 
larger energy range, (c) Phase space at the elliptical bound- 
ary with rotating (red lines from left to right) and oscillating 
(blue lines) regions and illustrations of trajectories. 
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the billiard boundary 8tt for which parallel trajectories 
show rotating motion. This is the case for the elliptical 
boundary of length C. Trajectories parallel to the hori- 
zontal boundary are precisely on the separatrix between 
oscillating and rotating motion. Therefore, integration 
over <5(sin 0) in eq. (|12[) gives half of the contribution for 
each of the two invariant regions of phase space and we 
have L Iot = C + r a /2. For the oscillating states we have 
L osc — ri, + r a /2. Note, that for the full ellipse the com- 
plete boundary belongs to the rotating region, L mt = 4£ 
and Lose = 0. Numerically we calculate the first 2568 
eigenstates of the desymmetrized elliptical billiard with 
r a = 1.0 and r& = 0.7 using the improved method of 
particular solutions 27]. They are characterized by the 
angular and the radial quantum number m and n. For 
each state we calculate the second constant of motion 
«mn [H|. If kJ,„ > f 2 the state is classified as rotat- 
ing and for n 2 nn < f 2 as oscillating. Figure EJa) shows 
the rotating and the oscillating spectral staircase for the 
elliptical billiard. We find excellent agreement with our 



prediction, eq. ([3]) (smooth solid lines). In fig. [3jb) we 
demonstrate that the boundary contribution of eq. ([3]) 
is in agreement with the difference of the numerical data 
and the first term of eq. ([3]) . We have confirmed that also 
under variation of the prediction eq. ([3]) agrees with 
numerics (not shown). 

A straightforward generalization of our results to Neu- 
mann boundary conditions is possible by changing the 
sign of the second term in eqs. ([2]), (|3]), (ITU1) . and (fT2l) . 
As interesting tasks there remains to find the higher or- 
der terms of dr(E) due to corners and curvature effects 
as well as to generalize our approach to systems with 
broken time-reversal symmetry and to three-dimensional 
cavities. Also the generalization of the approach to sys- 
tems with smooth potentials [32l ] is an open problem. 
Finally, it is now possible to study the spectral fluctua- 
tions around e?r (E) associated with a phase-space region 
r for generic billiards. 

We thank M. Sieber for valuable discussions and the 
DFG for support within the Forschergruppe 760 "Scat- 
tering Systems with Complex Dynamics" . 
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